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A crystal plasticity theory was developed for use in simulations of dynamic loading at high pres- 
sures and strain rates. At pressures of the order of the bulk modulus, compressions o(100%) may 
be induced. At strain rates o(10^)/s or higher, elastic strains may reach o(10%), which may change 
the orientation of the slip systems significantly with respect to the stress field. Elastic strain rather 
than stress was used in defining the local state, providing a more direct connection with electronic 
structure predictions and consistency with the treatment of compression in initial value problems in 
continuum dynamics. Plastic flow was treated through explicit slip systems, with flow on each sys- 
tem taken to occur by thermally-activated random jumps biased by the resolved stress. Compared 
with simple Arrhenius rates, the biased random jumps caused significant differences in plastic strain 
rate as a function of temperature and pressure, and provided a seamless transition to the ultimate 
theoretical strength of the material. The behavior of the theory was investigated for matter with 
approximate properties for Ta, demonstrating the importance of the high pressure, high strain rate 
contributions. 

PACS numbers: 62.20.F-, 62.20.D-, 62.50.Ef, 52.38.Mf 
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I. INTRODUCTION 



Plasticity has been investigated extensively in terms 
of slip systems for deformation at relatively slow rates, 
close to ambient pressure Some treatments incorpo- 
rate thermal activation of dislocation glide, but are usu- 
ally still constructed around an empirical yield stress for 
each slip system, and an empirical (power-law) harden- 
ing model [1, H, [j, [1| . In contrast, deformation at strain 
rates appropriate for shock loading is almost exclusively 
described in terms of isotropic elasticity and macro- 
scopic, empirical flow theories, usually treating plastic 
strain through a single, scalar parameter [tI. [sI. IqI. [lOj . 
These formulations are sub-Hookean, treating pressure 
and stress in a fundamentally different way which is an 
unnatural framework to extend to anisotropic elasticity. 
Plastic flow theories based explicitly on slip systems have 
only recently been developed in a form suitable for sim- 
ulating dynamic loading to moderate pressure (tens of 
percent of the bulk modulus), again using a flow rule for 
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each system based on a yield stress 

In principle, it is possible to predict any material prop- 
erties given an adequately detailed description of the 
composition and microstructure. At present, there are 
few properties that can be predicted as accurately as they 
can be measured by experiment: this is a deficiency in 
our understanding of condensed matter physics. Specifi- 
cally, it is not yet possible to predict the properties of 
polycrystalline 'engineering' materials to the accuracy 
where structures and devices can be constructed reliably 
and safely, without experiments to determine the mate- 
rial properties. One important application is the seed- 
ing of instabilities in inertial confinement fusion (ICF) 
implosions by spatial variations in the response of the 
thermonuclear fuel capsule In current ICF designs, 
plastic flow occurs during loading to several hundred gi- 
gapascals over time scales of a few nanoseconds: condi- 
tions that are far more extreme than considered in other 
applications involving plastic flow. Measurements of the 
velocity history and in-situ x-ray diffraction during shock 
loading of single crystals of Be on nanosecond time scales 
showed features more consistent with a pure thermal ac- 
tivation model rather than models based on a yield stress 
[T^ ]. so the present study focuses on the construction of 
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a crystal plasticity theory suitable for simulations of ex- 
treme conditions of pressure and strain rate using ther- 
mal activation rather than yield stresses. 

More generally, laser ablation is used increasingly in 
studies of the response of materials to dynamic loading, 
often exploring pressures to hundreds of gigapascals with 
durations of order 1 ns or shorter, and investigating the 
effect of crystal orientation and microstructure. Such ex- 
periments can only be interpreted in the context of high- 
pressure, high-rate crystal plasticity, and care is needed 
to consider the physics of deformation in adequate de- 
tail particularly when drawing conclusions or applying 
models on different (generally longer) time scales. 

It is reasonable to suppose that a quantitative un- 
derstanding of polycrystalline materials can be achieved 
through an understanding of the individual crystals mak- 
ing up the microstructure. The focus of the research re- 
ported here is the development of a theoretical approach 
to plasticity allowing the response of material microstruc- 
tures to be predicted given relevant properties of the sin- 
gle crystals, for deformation at high pressures and on 
short time scales. Materials with a grain size of order 
micrometers or less can be simulated for short periods 
using molecular dynamics (MD) [l4], but this is imprac- 
tical for the huge class of materials with larger grain sizes 
or where the response must be followed for long peri- 
ods. The interatomic potentials used in MD simulations 
are not always accurate, particularly for non-close-packed 
crystal structures. Dislocation dynamics (DD) simula- 
tions can treat larger volumes of material, but again the 
interaction and evolution rules for the dislocations are 
not necessarily accurate. The complementary approach 
adopted here is to simulate the response of polycrys- 
talline materials using continuum dynamics, where the 
microstructure is resolved and treated in terms of the 
response of individual crystals. As well as investigating 
physical processes occurring in dynamic deformation on 
length scales from single crystals to polycrystalline mi- 
crostructures, this type of model will help in the identi- 
fication of key material parameters and the development 
of physics-based models with reduced sets of parameters 
for continuum simulations in which the microstructure is 
not resolved. Resolved-microstructure simulations have 
been reported, and compared with experimental observa- 
tions [3, H, 01 , but not at strain rates and peak pressures 
high enough to apply to ICF. 

Because the conditions of interest are so far removed 
from the focus of previous research in crystal plasticity, it 
is important to preserve as much as possible of the sim- 
pler material models and numerical methods previously 
applied to such conditions; this constraint guided the 
structure of the crystal plasticity theory described here. 
We chose to retain the deviatoric decomposition of stress 
and strain used for isotropic strength models [15], allow- 
ing the isotropic pressure-compression relation to be cap- 
tured by a scalar equation of state (EOS), which is well- 
researched and reliable compared with the non-isotropic 
elastic and plastic response at high pressures and strain 



rates. Similarly, we followed an operator-splitting algo- 
rithm [isj for the simultaneous integration of elastic and 
plastic strain rates, widely used in the shock regime and 
suitable for stiff systems of equations. 

Pressures of the order of the bulk modulus induce 
compressions of tens of percent, which generally cause 
significant changes in elastic constants and plastic flow 
rates. Modern dynamic loading experiments commonly 
explore material response on nanosecond to picosecond 
time scales, where elastic strains may reach of order 10% 
p^ : under these conditions the orientation of slip sys- 
tems may change by several degrees even with no ro- 
tation of the material, so resolved shear stresses may 
change significantly. Here we present a theory of plas- 
tic flow derived consistently for high pressures and large 
elastic strains. Calculations of the importance of various 
contributions to the behavior of matter in these extreme 
conditions are demonstrated for Ta, which has been in- 
vestigated experimentally and using other theoretical ap- 
proaches, with many of the basic properties such as elas- 
tic constants predicted at high pressures. The purpose 
of this paper is not to present a complete and accurate 
theory for plasticity in Ta in the regime of interest, but to 
investigate novel properties of plastic flow under extreme 
conditions. 



II. CONTINUUM MECHANICS 

The problem of principal interest is to predict the re- 
sponse of spatially-resolved samples of known properties 
under the influence of a load which may vary with time 
and position. At each location in the sample at each 
instant of time, the material is described by its veloc- 
ity u and state s. This is an initial value problem: 
given {u, s}{f,to) over some region {r} S TZ'^ , what is 
{u,s\{f,t > to)7 Dynamic loading may be induced by 
the impact of material in different regions of the prob- 
lem, or by boundary conditions such as constrained val- 
ues of u{r,t) or the stress tensor T{r,t). for {r} S dTZ^ . 
(Note that, throughout this work, the stress is defined in 
the local rest frame of the material, i.e. it is the Cauchy 
stress.) In condensed matter on sub-microsecond time 
scales, heat conduction is often too slow to have a sig- 
nificant effect on the response of the material, and is ig- 
nored here. The equations of non-relativistic continuum 
dynamics describe the conservation of mass, momentum, 
and energy for matter fields with negligible diffusion of 
atoms, evolving under external and internal stresses. In 
Lagrangian form, i.e. along characteristics moving with 
the local material velocity u(r). 
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[15| where p is the mass density and i the specific internal 
energy. 

We follow the frame-based convention common in com- 
pressible hydrodynamics that 'Lagrangian' refers to a 
frame moving locally with the material and 'Eulerian' 
refers to a frame fixed in space. This is in contrast to 
the derivative-based convention common in continuum 
mechanics that 'Lagrangean' refers to derivatives with 
respect to the original, undeformed coordinate system 
and Eulerian' refers either to derivatives with respect 
to coordinates fixed in space or to coordinates moving 
locally with the material. Our meaning of 'Lagrangian' 
derivatives is with respect to coordinates translated along 
with the material, but not compressed or sheared. This 
convention leads to choices of tensor definitions which 
are somewhat different from those usually seen in contin- 
uum mechanics, but make it significantly easier to per- 
form integrations with respect to time over general large- 
deformation paths relevant to loading at high pressures 
and strain rates. 

The continuum dynamics equations are closed through 
the addition of constitutive relations describing the inher- 
ent properties of each material: {t(s); s(s, gradu, di/dt)} 
where s(r, t) is the field describing the local state of the 
material. The second law of thermodynamics provides a 
constraint on the constitutive model. r(s) is the EOS; 
s{s, gradu, di/dt) describes the evolution of the state un- 
der an applied deformation and heating rate. As well as 
experiencing compression and work from mechanical de- 
formation, s(r, t) can evolve through internal processes 
such as plastic flow. Subsuming Dp/Dt and Di/Dt into 
the evolution relation, 

^^^ = s[s{r,t),U{r,t)] : [/ = grad u(f, t) . (4) 

The crystal plasticity model described here was imple- 
mented in a software package of material models, pro- 
viding a common interface for different varieties of con- 
tinuum dynamics program [TtI . [isj . The software struc- 
ture is object-oriented, the continuum dynamics program 
storing the local material state s regardless of type, and 
the material package performing calculations including 
t(s) and the rate of change of state given grad u or di/dt. 

There is an inconsistency in the standard contin- 
uum dynamics treatment of scalar (pressure) and tensor 
(stress) response. The scalar EOS expresses the pressure 
p{p, T) as the dependent quantity, which is the most con- 
venient form for use in the continuum equations. Stan- 
dard practice is to use sub-Hookean elasticity [llj, in 
which the stress tensor is found by integration 

T — ce (5) 

where c are elastic constants and e is the elastic strain 
rate tensor. Thus the isotropic and deviatoric contri- 
butions to stress are not treated in an equivalent way: 
the pressure is calculated from a local state involving a 
strain- like parameter (mass density), whereas the devia- 
toric stress is part of the local state and evolves by time- 
integration. This is an undesirable distinction which 



complicates the problem. The stress can be considered 
as a direct response of the material to the instantaneous 
state of elastic strain: T{e,T). This relation can be pre- 
dicted directly with electronic structure calculations of 
the stress tensor in a solid for a given compression and 
elastic strain state [l^, and is a direct generalization of 
the scalar EOS. 

The constitutive response of a crystal is described most 
naturally in terms of a coordinate system fixed with re- 
spect to the symmetry directions of the crystal, such as 
its lattice vectors. The orientation of elements of mate- 
rial varies with position and time in the problem, and is 
described completely by a local rotation between crys- 
tal and problem coordinates. Thus the state s is chosen 
to give the elastic strain e in the coordinate frame of 
the crystal, and includes the rotation matrix R from the 
crystal to the problem coordinates. This is a polar de- 
composition of the instantaneous deformation gradient 
F = R^V, where e is related to the right stretch ten- 
sor V through a specific choice of finite strain measure. 
The evolution equations for V (hence e) and R implic- 
itly include contributions from the 'external' evolution of 
the continuum and from the 'internal' evolution of the 
local material state by plastic flow: subscripts c and p 
respectively: 

V = V^ + Vp^ e = ec + ep, R = Rc + Rp. (6) 

Integrating over time, this is a multiplicative decomposi- 
tion of the deformation gradient F into elastic and plastic 
contributions, though only the elastic contribution is re- 
quired explicitly. 

At any instant of time, the continuum dynamics cal- 
culation provides gradu(r, i) = U. This is decomposed 
into symmetric and antisymmetric terms S and A, re- 
spectively describing strain and rotation rates: 

U = S + A : S =^{U + U^),A=^{U -U'^). {7) 

Neglecting plastic flow for the moment, the rate of change 
of the right stretch tensor V from the evolution of the 
continuum is 

V = RSR^V. (8) 

The external rotation rate is found similarly by noting 
that the incremental rotation caused by A over a short 
increment of time is 5Rc ~ I + 5tA^ to order 5t, so 

Rc = A^R + RA. (9) 

The strain tensor e is calculated from V , or Cc from Vc, 
through the use of a finite strain measure. The main con- 
sideration in the choice of strain measure is consistency 
with the measure chosen when inferring the stress tensor 
T from measured or calculated states of finite strain. Pos- 
sible choices include In V and {V"^ — I)/rn, where m = 2 
gives the Green-St Venant measure [20|. The strain is 
objective, which makes the formulation computationally 
convenient. 
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The formulation used for the constitutive model was 
that, for a given composition and crystal structure of a 
material, in the frame of reference of the crystal, which is 
rotated with respect to the continuum dynamics frame, 
the stress tensor depends only on the elastic strain 
tensor e and the temperature T. Strictly, also de- 
pends on the defect structure in the material; for most 
situations of practical or research interest, the vast ma- 
jority of atoms in a crystal are in an environment which 
has the underlying symmetry of the crystal, so the direct 
contribution of defects to was ignored. Plastic flow 
acts to reduce shear stresses, i.e. to reduce the shear 
components of elastic strain. Plastic flow is mediated 
by the motion of defects such as dislocations, defined 
with respect to the instantaneous crystal axes. At the 
level of the crystal lattice, the plastic state comprises the 
population of defects {<j>} on each slip system and the 
population of obstacles to defect motion {w}. The right 
stretch tensor is needed in order to relate the defect lat- 
tice vectors to the instantaneous elastic strain state. For 
the strain measures considered, it is more efficient to cal- 
culate e{V) than V{e), so V was used in defining the 
material state. Thus the local, instantaneous state of 
the material is s{f,t) = {R,V,T, {(j)}, {oj}}. The consti- 
tutive relation describes Tx{s), the relaxation of elastic 
strain by plastic fiow ep(s) and the resulting conversion 
of elastic to thermal energy T{s, Cp), and the evolution of 
the plastic state as flow occurs: {0}(s, Cp) and {w}(s, Bp). 
A hyperelastic stress-strain relation was used: 

Tx = c{s)e (10) 

where c(s) is the tensor of elastic 'constants,' explicit 
functions of elastic strain (particularly compression) and 
temperature. It is easier to ensure that the stress is a 
path-independent function of elastic strain with a hyper- 
elastic relation than a hypoelastic relation 

Ta, = S(s)e. (11) 

The (objective) stress in the frame of the problem is 
found using the local rotation: 

T(f, t) = R^{f, t)T.,{f, t)R{f, t). (12) 

The effect of plastic flow can be calculated and inte- 
grated simultaneously with the flow equations through 
s(s, gradu, 9«/9t). For time- integration by finite differ- 
ence numerical schemes, the natural time increment from 
the continuum dynamics may be very different from that 
required for plastic flow. Similarly, plastic heating may 
change the state enough that numerical schemes integrat- 
ing forward in time become unstable. For this reason, for 
most purposes plastic flow was operator-split from the 
continuum dynamics (isj : the continuum equations were 
integrated at constant plastic state, then plastic flow was 
calculated at constant total strain. For a complicated 
material response in which the rate of change of general- 
ized state is a sum of logically-distinct contributions 

s = ^Sj, (13) 



operator-splitting is the use of the approximate integral 
relation 

/t + At nt + At 

s{s)dtc^Yl s^{s^)dt (14) 

where Si is the result of the previous integrations j < i, 
i.e. a sequence of integrations over the full time increment 
At in which all but one of the contributions in turn are 
ignored. This decomposition is used widely in continuum 
dynamics involving shocks and other high-rate flows, for 
example to decouple Eulerian remaps in different dimen- 
sions from each other and from Lagrangian integrations 
, and to decouple stiff chemistry from the continuum 
in reactive flow [23|. It is typically second-order accu- 
rate in time. A previously-reported algorithm employed 
a plastic predictor and an elastic corrector [22] and was 
demonstrated to be stable and accurate. However, this 
is not a natural framework for plastic flow as a response 
to dynamic loading at extremely high rates and compres- 
sions, as the material response is best regarded as elastic 
with plasticity occurring as a stress-relaxing mechanism, 
as demonstrated by observations such as the deca y o f 
elastic precursor waves in shock loading [2^ [U, [2^ 12^ . 
Our formulation also appears to be more suitable for de- 
viatoric strength with an EOS in the strong shock regime. 
Operator-splitting allows the plasticity algorithm to be 
more stable than is generally possible when integrated si- 
multaneously with explicit-time hydrodynamics (because 
plastic flow cannot then lead to time increments that are 
inconsistent between the continuum and the local plastic 
state) with little change in accuracy and similar or better 
computational efficiency. When integrating the contin- 
uum equations with operator-splitting, the constitutive 
model was split into a series of partial integrations 

Se = s(t) + / Se{s,gTa,du,die/dt) dt (15) 

Jt;s{t) 

Sh = Se+ Sh{s,0,dih/dt)dt (16) 

Jt;s<, 

{t')+ / Spis,0,dip/dt)dt, (17) 

where At = t' — t is the finite time increment, Se is the 
state evolution under the infiuence of the continuum ve- 
locity gradient with no plastic flow and heating die/dt 
restricted to isentropic compression plus kinematic vis- 
cous forces, Sh is heating dih/dt from external sources 
(including radiation deposition, conduction, and artificial 
viscosity used to stabilize shock waves) , and Sp is plastic 
relaxation with the total strain held constant with heat- 
ing dip/dt only from dissipation of relieved elastic strain 
energy. Integration proceeds by taking s(t') = Sp. 

The software implementation made wide use of object- 
oriented techniques, particularly polymorphism and in- 
heritance. This allows the constituents of the material 
models to be chosen when the program is run, in partic- 
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ular the form and type of the EOS, and the type of gen- 
eral functions such as polynomials or tabulations. The 
ability to define functions in a general way in the input 
to a program, rather than having to be specific in the 
source code, is convenient for representing the evolution 
of defect populations in the crystal (including disloca- 
tions and pinning sites) as it is less clear which precise 
fmictional forms to use. 



III. DEVIATORIC ELASTICITY 

In high pressure loading, plastic fiow is usually fast 
enough that the isotropic pressure has greater magnitude 
than shear components of stress. As discussed above, 
unusually large elastic shear strains may reach ~10%, 
and are more commonly At pressures the order of 

100 GPa, isotropic compressions are typically several tens 
of percent, so the isotropic stress clearly dominates. For 
this reason, we use the deviatoric formulation of elastic- 
ity. The stress and strain tensors r and e are decomposed 
into isotropic and deviatoric contributions: 

e{f,t) = €{f,t) + ^i{r,t)I : ^^^Tre (18) 

o 

t(s) = a(s)-p(s)/ : p=-iTrT. (19) 

The mean pressure p{s) is represented by a scalar EOS, 
in practice using the mass density p{f, t) rather than the 
compression nif^t) in the state s{r,t). The deviatoric 
stress-strain relation is 



a{s) = C{s)e 



(20) 



[27| where C is the tensor of deviatoric elastic 'constants,' 
explicit functions of p, T, and in principle e. C are ob- 
tained from c by subtracting the iscntropic bulk modulus 
B from all elements affecting the normal components of 
stress. Here we take 



B = 



dp 
dp 



(21) 



which is well-defined for all crystal systems. (In general 
for non-cubic systems, e 7^ for cr = 0.) Voigt notation 
was used for strain and stress tensors, so c and C were 
6x6 matrices. In Voigt notation. 
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(22) 



Non-linear elasticity is accommodated by allowing C or 
c to be a function of the elastic strain deviator - this is 
potentially important for elastic strains above the percent 
level. 



Strictly, the scalar EOS is also a function of the strain 
deviator [16|. Ignoring this dependence may lead to 
0(1%) errors at the elastic strains considered here. This 
refinement will be considered further in future work. 

Thermodynamic consistency depends on the structure 
of the constitutive model. In thermoelasticity, the en- 
tropy at any strain is a function of temperature only. 
Once the defect structure associated with plasticity is in- 
cluded, it also contributes to the entropy. However, de- 
fects generally affect only a small fraction of the atoms. 
Therefore, above the Debye temperature of the material, 
the defect entropy is small compared with the lattice- 
vibrational entropy. In this work, we neglected the en- 
tropy of the defects, and used thermodynamically consis- 
tent scalar EOS to capture the relation between specific 
internal energy i and temperature T. For consistency 
with the construction of the scalar EOS, i was chosen to 
exclude the elastic distortional energy, which was stored 
and conserved implicitly through the elastic strain. 

The scalar EOS has been obtained for many materials 
[HjII^,!!^!, by high-pressure experiments such as impact- 
induced shock waves fsi'l, and from electronic structure 
calculations J2]. Elastic constants have been measured 
at STP for many materials, and can also be predicted as 
a function of compression, temperature, and strain from 
electronic structure calculations [19l|. We took one such 
prediction for Ta up to 1 TPa [33| . The reported varia- 
tion of (Cii — Ci2)/2 with pressure was used to calculate 
cii(p) and Ci2(p), using a semi-empirical scalar EOS [2§| 
to obtain pressure p(p) and bulk sound speed c^(p) along 
the zero-kelvin isotherm. The elastic constants predicted 
at zero pressure deviated by around 10% from experi- 
mental measurements [33|. 

Conventionally, elastic constants C are expressed as 
functions of pressure p and temperature T [y, |Tl[. In 
normal materials, dC/dp\T > and dC/dT\p < 0, with 
\dC/dp\T\ > \dC/dT\p\. It is more convenient to calcu- 
late C{p,T) from electronic structure calculations rather 
than C{p,T), and generally \dC/dT\p\ < \dC/dT\p\. 
In fact, in some cases it may be reasonable to ignore 
the temperature-dependence of C{p,T). For instance, 
considering the variation of the polycrystalline shear 
modulus G for Ta [30], dG/dp\T = 1.4510-2Go/GPa 
and dG/dT\p = -I.SIO-^Gq/K, implying dG/dT\p = 
— 1.29810~^Go/K. The decomposition in terms of p and 
T is therefore at least as appropriate for elasticity over a 
wide range of compressions. 



IV. CRYSTAL PLASTICITY 

At the level of a single crystal or grain, plastic flow 
is mediated by defects - dislocations or disclinations - 
which act along a set of slip systems defined with re- 
spect to the lattice vectors. The theory described here is 
based on previous work on rate-dependent crystal plas- 
ticity [H , but generalized to treat high compressions and 
elastic strains, which entails the inclusion of some cor- 
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rections and refinements. The theory was developed for 
dislocation-mediated plasticity, but is expressed in a form 
which enables it to be used to represent flow by twin- 
ning, in at least a simplified manner. Where appropriate, 
we refer to defects rather than dislocations, to make the 
point that much of the theory should apply to disclina- 
tions. 

In general, twinning and slip are potentially compet- 
ing mechanisms which interact with each other. The ori- 
entation of dislocations and slip systems changes when 
material twins. The presence of dislocations can impede 
the motion of a twin boundary, the structure and mobil- 
ity of a dislocation are altered if it passes through a twin 
boundary, and the presence of twinned material alters 
the strain field from, and forces on, a nearby dislocation. 
The treatment of plasticity adopted here is to resolve 
each element of material with a single crystal orienta- 
tion, populated by defects described by a set of homoge- 
neous densities, mediating strain by their motion. When 
applied to twinning, this approach lends itself to sim- 
ulations where either the evolution of twins is resolved 
spatially, or strain occurs only in a single plane, or the 
volume fraction of twinned material is small. These limi- 
tations could be removed by allowing volume fractions of 
each variant of twin orientation, as has been reported pre- 
viously 35, 36, 37; 38, 39]. For practical reasons of rapid 
application to real problems, previous models of twinning 
have been constructed around a yield stress for each de- 
formation mechanism, with modifications to power-law 
hardening, rather than explicit defect kinetics. Simplifi- 
cations have been made such as assuming that dislocation 
motion would be negligible in twinned regions. Gener- 
ally, models of twinning are less advanced than those of 
crystal plasticity, and we present the dislocation-focused 
work here as a starting point for future developments in 
twinning rather than a complete framework for a model. 

The local state in the material is described by a set of 
defect populations and a set of obstacles to defect mo- 
tion. Each defect population is characterized by its type 
(e.g. edge, screw), orientation with respect to the crystal 
lattice vectors (e.g. by its Burgers vector b), and local 
density (j). In reality, extended defects (dislocations and 
disclinations) are curved, and may be of mixed type. In 
this work, mixed-type or mixed-orientation defects are 
accounted for by assigning them proportionately to the 
appropriate 'pure' populations. The local state of the 
material thus comprises the sets of defect and obstacle 
densities, {</>&} and {w}. For each population of defects, 
there is a set of planes and directions {n, d} in which the 
crystal shears. Each specific combination of 6, n, and d 
is referred to as a slip system, but different slip systems 
share the same defect population if b is the same. There 
are separate defect populations with 6 and —6, but for 
given b there are not separate slip systems with direction 
n and —n or d and —d: the latter is treated as reverse 
motion on the slip system with direction d. The obsta- 
cles experienced by defects with b are the same as those 



with —b, so the local state of the material was taken more 
compactly as 0^}; {w}. 

At the atomic level, the relevant densities of defects 
and obstacles are expressed with respect to the densities 
of atoms. Thus the defect and obstacle densities do not 
change with elastic compression. 

Under the action of the local stress, the defects move 
thus causing plastic flow which relieves elastic strain, and 
the populations of defects and obstacles evolve. The plas- 
tic strain rate was calculated for each slip system inde- 
pendently. The time increment for integrating plastic 
flow was chosen based on the total plastic strain rate, 
so that the increment of plastic strain was limited to a 
small value (such as 0.01) per time step. Since plastic 
flow was operator-split from the continuum dynamics, it 
was thus subcycled at constant total strain over the time 
step required by the continuum dynamics. This approach 
is used widely for the simultaneous solution of stiff equa- 
tions with continuum dynamics, for example in reactive 
flow [2l| . In this instance, it greatly simplifies the calcu- 
lation of potentially competing slip systems. 

Given the stress tensor t, the stress resolved along the 
slip system is 



— d ■ a ■ n = Ur 



(23) 



the resolved deviatoric stress, since n and d are always 
orthogonal. In the frame of the undistorted crystal, slip 
in a given system produces a change of total strain of 



Act oc —d ® n, 



(24) 



which is traceless - no volume change - since n and d 
are orthogonal. The strain increment can be split into 
elastic (unloading) and rotational parts, 



Acp — Aep 



-i {Act 



Ael) 



Ai?„ 



'2 i^^T 



Ae^) 



(25) 
(26) 



where the sign of ARp is opposite to that of Aep because 
R rotates from the crystal coordinates to the instanta- 
neous orientation in the laboratory space. The conven- 
tion chosen here is to calculate the strain increment with 
respect to the undistorted lattice, and accumulate spe- 
cific or per-atom defect and obstacle densities. If Hq and 
d'g are the normal and direction as real space vectors with 
respect to the undeformed crystal, in the deformed crys- 
tal 



Vd'o, 



(27) 



where the latter equation uses Nanson's formula. The 
change in elastic strain is then calculated with respect to 
the distorted lattice, using the right stretch tensor V, 



Act ^ d' i^n' 



(28) 
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and used to calculate plastic flow increments of stretch 



AV^ and rotation Ai?p: 



AVp = --(ASt + A. 
ARp = +i (Aer - Ae"^ 



(29) 
(30) 



This correction takes into account the orientation effect 
of finite strain and the scaling of the lattice vectors by 
isotropic compression (Fig. [T|). 

In Ta under high-rate loading conditions, screw dis- 
locations dominate the behavior [s^. There are four 
distinct dislocation populations (and two signs of each) , 
with 6 = i(III). Each population with ±6 may slip along 
{110} planes (3 systems), {112} planes (3 systems), and 
{123} planes (6 systems) (see e.g. [i^ for a full enumer- 
ation of the Burgers vectors and planes). For a crystal 
compressed uniaxially along [100], the effect of finite elas- 
tic strain is an o(10%) effect for strains of tens of percent 
(Fig.©. 

Defect motion proceeds as a thermally-activated hop 
past the Peierls barrier. In principle, hopping may occur 
in the direction of increasing shear energy; it is more 
likely to occur in the direction of decreasing shear energy. 
With no elastic stress, hopping is equally likely in either 
direction. An applied stress biases the average hopping 
direction by altering the effective Peierls barrier by the 
work that would be done against the applied stress at the 
peak of the barrier. If the Peierls barrier energy is E},, the 
effective barrier is E{, ± AE where AE = arfl/2, I is the 
interatomic spacing in the direction of slip (a function 
of compression), and / scales arb to the units used to 
express Ei,, e.g. an area per atom if Ei, is expressed as 
the energy per atom. The net hopping rate A should 
thus follow an Arrhenius form with competing forward 
and reverse contributions, 



X = Z{p,T) 



Ei,+AE 



(31) 



where Z is the rate at which the defect attempts to hop 
- the vibration frequency of the atoms - and fcs is Boltz- 
mann's constant. Rearranging to factorizc out the un- 
strained Arrhenius barrier for single-direction hopping. 



A = 2Z{p, T) exp 



^^%inh^. (32) 



Each exponential term is the probability of a successful 
hop: the Arrhenius form applies when the probability is 
less than one, and each exponential was limited to one. 
Strictly, the Arrhenius form applies when the hopping 
attempts are described by a Boltzmann distribution. In 
a solid, the distribution is given by the population of 
phonon modes, which gives a modified probability par- 
ticularly below the Debye temperature, when the popu- 
lation of the phonon modes falls significantly and zero- 
point energy becomes more evident. The hyperbolic sine 



dependence has been discussed previously in the context 
of creep 41], though the treatment of applied strain and 
stress is different here. When dislocations move, the rate 
at which a straight dislocation develops a forward-moving 
kink is lower than the rate at which the kink spreads lat- 
erally as adjacent sections of the dislocation move for- 
ward, because the energy of an isolated kink is higher. 
On average, the net successful hopping rate may depend 
on hops of the same direction by several adjacent atoms 
{N say), so 



\^Z{p,T) 



-N- 



-N- 



(33) 



If there were no obstacles, the overall plastic strain rate 
on the slip system would be 



a.bnd = hnd {4>b + <l>b) ■ 



(34) 



The effect of obstacles is to reduce the mobility of dislo- 
cations and hence reduce the mean plastic strain rate, be- 
cause of the finite energy required to generate additional 
lengths of curved dislocations. This energy includes con- 
tributions from the dislocation core and the strain field, 
which for a curved dislocation generally includes self- 
interactions. Obstacles include point defects (vacancies, 
substitutional impurities, and interstitials), other popu- 
lations of dislocations, and disclinations. If there were 
only one type of obstacle (density lo) and it were able to 
pin dislocations perfectly, the resulting strain rate on the 
slip system would be 



(35) 



where P{Q) ~ 1 and P(l) — 0. Here we choose P{lo) — 
1 — uj for < < 1 per atom. Each type of obstacle has 
a different inherent effectiveness in reducing the mobility 
of dislocations of each type, represented here by a factor 
f^ibnd for the zth type of obstacle. Dislocations may be 
able to pass obstacles, which can be characterized by an 
Arrhenius rate so that the overall effectiveness of a given 
type of obstacle is 



hnd 



13. 



hnd 



1 



(36) 



where Eo is the energy barrier for the dislocation to pass 
the obstacle and AE is calculated as for Peierls barriers. 
Here we assume that the density of obstacles is much 
lower than that of the atoms in the crystal, so reverse 
hopping should be neglected though barrier lowering is 
included. Combining the effects of all obstacles, 



abnd = Afcnd (06 -I- 4>l) P{uJbnd), 



where 



l^bnd 



'ibnd^i- 



(37) 



(38) 



Each slip system contributes to the unloading of the elas- 
tic strain at a rate in each system of 



f-bnd 



1 . / 

--jOibnd y 



d ® n + n ® d 



(39) 




FIG. 2: Effect of finite elastic strain on shear stress resolved 
along slip systems in Ta, for uniaxial compression along [100] . 
Reference value is resolved shear stress on each system ne- 
glecting the effect of finite strain. 

As plastic flow is expressed as a strain rate which is a 
function of the instantaneous, local strain state, this for- 
mulation can be regarded as viscoplastic. 

In previous work, the attempt frequency Z for the 
motion of a dislocation has been held to be much less 
than the attempt frequency for an atom (ssj . How- 
ever, vibrations of the atoms provide the fundamental 
mechanism for dislocations to move. We note that the 
reverse-hopping contribution and the multiple-atom fac- 
tor decrease the net successful hop rate, and could be 
interpreted as a smaller attempt frequency in an over- 
simplified Arrhenius treatment. 

In this theory, the effect of the applied stress is to tilt 
the energy surface seen by atoms constituting the defect 
as they hop in the slip direction. The ultimate strength of 
the material is given by the maximum gradient in the un- 
stressed energy surface: when the resolved stress reaches 



FIG. 3: Effect of resolved stress on the periodic potential 
experienced by a moving defect (displacement in units of the 
interatomic spacing, amplitude in units of the Peierls barrier, 
stress in units of the ultimate strength). An applied stress 
reduces the barrier height in the direction of (mean) plastic 
strain; at the ultimate strength the barrier disappears in that 
direction. 

this level, there is no barrier to defect motion. For a sine 
wave potential, 

y{x) = —sm—, (40) 

the ultimate strength is nEb/l. (Fig.[3l) 

The reverse-hop term has not been included in other 
plasticity studies using Arrhenius-based flow rates. Its 
effect is greatest at low strains, i.e. at low strain rates, 
where the stress biasing term is smallest. For example, in 
Ta (flow stress observed in microsecond-scale loading ex- 
periments '^l GPa [lOl, distance between adjacent Peierls 
valleys I ~ 1.5a where a = 3.308A), the effect of the 
reverse term is to reduce the strain rate by tens of per- 
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temperature (K) temperature (K) 



FIG. 4: Effect of the barrier-lowering and reverse-hopping 
terms on the strain rate, for Ta at zero pressure with 1 GPa 
resolved stress. 



FIG. 5; Predicted effect of resolved shear stress (GPa) and 
temperature on the average speed of (111){110} dislocation, 
for Ta at 16.7 g/cm^ 



cent for temperatures over a few hundred Kelvin, and 
the lowering of the barrier by the applied stress increases 
the strain rate by several tens of percent. Comparing the 
overall effect with the nominal Arrhenius hop rate, the 
net effect is to increase the strain rate by hundreds of per- 
cent at temperatures below about 310 K and to decrease 
it by tens of percent at higher temperatures. (Fig. [4]). 

The effect of stress and temperature can also be seen 
by calculating the average speed of a dislocation moving 
through the lattice, 

Ud = A|6| (41) 

For temperatures of a few hundred kelvins and resolved 
shear stress (rss) well below the Peierls stress, the disloca- 
tion speed increases with temperature and rss, and has a 
value of similar order to kinetic Monte-Carlo predictions 
using an explicit representation of the double-kink mech- 
anism [i^]. However, once the rss exceeds the Peierls 
stress, the dislocation speed is limited to a value around 
half of the bulk sound speed, because of saturation of the 
probability for forward jumps. At high temperatures, the 
average dislocation speed decreases with temperature as 
the probabihty for reverse jumps increases. (Fig. [3) 

Repeating these predictions for Ta compressed to 
20 g/cm^^, the same trends are evident, but the saturation 
stress is greater than 6 GPa and the maximum disloca- 
tion speed is higher in proportion to the rise in frequency 
of atomic vibrations (Fig. [6|). 

The total working rate \s W = \\<jip\\ where Cp = 
— J2bnd^bnd- This is the rate at which elastic strain 
energy is relieved; it is not all converted into heat as 
some remains stored as the energy of the defect struc- 
tures. The propor tion fp converted to heat is thought to 
be 85-95% [43, |4J|; this was treated here as a parameter 
in the theory, but could be calculated by subtracting the 



3000 




temperature (K) 

FIG. 6: Predicted effect of resolved shear stress (GPa) and 
temperature on the average speed of (111){110} dislocation, 
for Ta isentropically compressed to 20g/cm'* (around 49 GPa 
on the principal isentrope). 

core and interaction energy of the defect population, say 
E = ^r]b(f)b + ^Vbb'(l)b<Pb'- (42) 

b bb' 

Such a calculation has been attempted previously, us- 
ing a yield stress and power-law hardening model for the 
crystal, but so far without an accurate treatment of the 
scalar EOS or elasticity at high compressions 

Since the rate of plastic flow is as a sum over slip sys- 
tems, and the rate is calculated independently for each 
system, this theory automatically treats non-Schmid be- 
havior where flow may occur on a slip system which does 
not have the greatest rss ^Bl, if it is kinetically favorable. 
This condition may occur if the system with greatest rss 
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has insufficient dislocations to relieve the elastic stress, 
or if the Peierls barrier is higher than for a competing 
system of lower rss. As the rate of plastic flow on a slip 
system is non-linear in rss, when competing slip systems 
have the same defect populations and kinetics, flow is 
dominated by the system with greatest rss, giving nor- 
mal Schmid behavior. 

As plastic flow occurs, the populations of dislocations 
and obstacles evolve. Sections of dislocation may anni- 
hilate if they reach a grain boundary or a dislocation of 
opposite Bur gers vector. In most situations studied ex- 
perimentally [4^, the dominant effect is the creation of 
additional lengths of dislocation by Frank-Read sources 
between obstacles [47[- Dislocations may be pinned by 
intersection with other defects in the lattice such as dis- 
locations in other slip systems. In the present work, the 
dislocation density was taken to evolve according to the 
generic relation 



[Fi{(l)b,ujb) - F2{(l)b<Pb)] ctb 



where 



E 

n,d 



(43) 



(44) 



the Fi are respectively spontaneous dislocation genera- 
tion, net dislocation generation by flow (including an- 
nihilation with the grain boundary), and annihilation 
with dislocations of opposite Burgers vector, and the 
Rij are reaction terms between dislocations of different 
types, which may involve Arrhenius rates to cross bar- 
riers. Spontaneous generation of defects is a reversible 
process limited by an energy barrier Eg with a rate 



Fo = Z{p,T) 



(45) 



Compared with flow by defect motion, spontaneous gen- 
eration can affect any atom in the crystal, but has a high 
energy barrier. 

Dislocation reactions may make a significant difference 
to the plastic flow rate. Mobile dislocations can combine 
to form sessile dislocations, reducing the density of mo- 
bile dislocations and acting as obstacles to the remaining 
mobile dislocations. 

Obstacles which are dislocations on other slip systems 
are described by the dislocation density on those slip sys- 
tems. In this case, the effectiveness parameters Pibnd 
constitute a hardening matrix. 

The vibration rate Z{p, T) can be estimated from elec- 
tronic structure calculations or measured for instance by 
low-energy neutron scattering. The Peierls barrier can 
also be estimated from electronic structure calculations 
with appropriate variations in elastic strain, or inferred 
from interatomic potentials. The initial density of de- 
fects and obstacles or pins is best obtained by analysis of 
the microstructure. The functions describing the evolu- 
tion of dislocations and obstacles, including the harden- 
ing parameters Pibndi can be estimated from MD or DD 



simulations. The fraction of plastic work converted to 
heat could be estimated from these simulations. 



V. EXAMPLE SIMULATIONS OF PLASTIC 
FLOW IN TANTALUM 

As an illustration of the behavior of constitutive mod- 
els based on the theory described above, we consider the 
time-dependent relaxation of shear stress in material sub- 
jected to an instantaneous uniaxial compression. This 
scenario represents the conditions induced by rapid load- 
ing such as laser ablation on a single crystal. It has been 
observed experimentally in a variety of materials includ- 
ing Be [l3|, and is related to the phenomenon of decay 
of a supported elastic wave [i^. As before, the material 
was taken to be 'Ta-like,' using parameters taken from 
the literature. 

The elastic constants and slip systems were as de- 
scribed above. The frequency factor Z(p, T) for defect 
motion was estimated from the ambient Debye temper- 
ature of Ta. Z{p) was assumed to be dominated by the 
variation of the elastic constants, so Z ^ B/ p, i.e. 
proportional to the bulk sound speed. Any temperature 
variation of Z was neglected. The Peierls barrier Eb{p^ T) 
to defect motion was taken from the zero-pressure value 
of O.lSeV calculated for rigid shear without lattice re- 
laxation [3^ and the calculated variation of the ultimate 
strength Tc of Ta with pressure. 



^.0.12^ 
dp dp 



(46) 



|3a | where G is the shear modulus. dG/dp was estimated 
from the pressure-hardening term in the Steinberg- 
Guinan strength model i^SOl] , assuming 



1 dEb 

Eb dp 



1 cLtc 
Tc dp 



(47) 



As with the frequency factor, any explicit temperature- 
dependence of the Peierls barrier was neglected. In both 
cases, there was an implicit temperature dependence as 
functions of p, as they were expressed in terms of p. 

The initial dislocation density in any crystalline ma- 
terial depends completely on the details of preparation: 
cooling rate from a melt, amount of subsequent work- 
ing, temperature at which working occurred, and so on. 
For metal crystals grown from molten material, e^^ by 
zone melting, the dislocation density is typically [49| in 
the range 1 to 10 per /xm^, which is equivalent to 6 to 
60 X 10~^ per atom. The population of each orienta- 
tion depends on, for example, the temperature gradient 
with respect to the lattice vectors during the formation 
of the crystal, and the orientation of the principal axes 
of the strain fields applied during working, again with re- 
spect to the lattice vectors. The dislocation density was 
taken to be distributed uniformly over all orientations of 
dislocation. Obstacles include dislocations on other slip 
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FIG. 7: Simulated stress relaxation in a Ta crystal compressed 
elastically by 1% along [100] , with no obstacles and no evolu- 
tion of the dislocation densities. The dislocation density was 
IQ/ jsre? (solid) and 100/pim'^ (dashed), distributed uniformly 
over all systems. 



systems, grain boundaries, impurities, and interstitials, 
and also depend completely on the details of manufac- 
ture. Simulations were performed with different obstacle 
densities, and different models of obstacle evolution. The 
hardening parameters Pibnd followed previous studies [1] 
in taking /3 = 1 for self-hardening and (3 ^ \ A for cross- 
hardening. Dislocation reactions were not considered. 

Simulations were performed of Ta crystals compressed 
along [100] and [110]. We consider two general classes of 
simulation, which explore different aspects of the mate- 
rial response: an initial elastic strain followed by plastic 
relaxation (which closely mirrors the characteristics of 
shock- loading experiments), and deformation at a con- 
stant strain rate (which is equivalent to more conven- 
tional dynamic loading tests). Except for the study 
of initial temperature, the sample was initially at STP. 
The orientation was chosen so that compression occurred 
along the x-direction in the coordinate reference frame. 



A. Elastic strain followed by plastic relcLxation 

Following the initial, instantaneous compression, the 
normal and transverse components of stress, th and 
T22 ~ 733 respectively, relaxed toward the mean pres- 
sure (Fig. [7]). As one would expect from the structure of 
the model, increasing the obstacle density decreased the 
rate of stress relaxation, and the evolution (increase) of 
the defect density by interaction with obstacles increased 
the rate of stress relaxation. However, the detailed choice 
of the functional forms and parameters, which is mani- 
fested as work-hardening, is speculative without guidance 
from MD and DD simulations, and will be described sep- 
arately. 
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FIG. 8: Effect of initial compression on simulated stress re- 
laxation in a Ta crystal. The initial compression was along 
[100], by fractional amounts shown in the graph. The dislo- 
cation density was held constant at IQ/ fxm^ , and there were 
no obstacles. Stress components rn (upper) and T22 (lower) 
are shown for each initial compression. 



Simulations were performed of plastic relaxation for 
different initial elastic strains. For initial strains greater 
than 5%, the barrier-lowering effect of the applied stress 
gave an abrupt decrease in the time required to relax 
to a state of isotropic stress (Fig. [8]). Simulations were 
performed for different values of the initial temperature 
from 100 to 1000 K. For simplicity, the temperature was 
applied at the same starting density, which gave a varying 
mean pressure. The time to relax to an isotropic stress 
state showed a clear dependence on temperature (Fig. [9]) . 

The stress relaxation varied significantly between [100] 
and [110] orientations parallel with the loading direction. 
For loading along [110], the initial stress was lower and 
relaxation took longer. There were also significant differ- 
ences in the other stress components (Fig. [TU]) . 



B. Constant strain rate 

During deformation at a constant strain rate, the finite 
rate at which defect-mediated plasticity could relax shear 
stresses allowed a finite elastic stress to be maintained, 
i.e. a finite fiow stress. Simulations were performed for 
deformation at different strain rates, held constant for 
a period of time. Strain rates were applied for uniaxial 
compression in the x-direction (representing ramp waves) 
and for pure shear in the xy plane. In uniaxial com- 
pression, the shear stress increased with compression at 
high strain rates and strains because the elastic constants 
increased with compression. The shear stress exhibited 
oscillations caused by competition between thermal soft- 
ening from plastic heating and pressure hardening from 
compression. For shear deformation, the flow stress de- 
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FIG. 9: Effect of initial temperature (increments of 100 K) on 
simulated stress relaxation in a Ta crystal. Heating was iso- 
choric, so the different initial temperature gave different initial 
pressures. The initial compression was 1% along [100]. The 
dislocation density was held constant at 10//.im^ and there 
were no obstacles. 
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FIG. 10: Effect of orientation on simulated stress relaxation 
in a Ta crystal. The initial compression was 1% along the 
X-axis. The dislocation density was held constant at 10//im^ 
and there were no obstacles. All non-zero stress components 
are shown for a crystal with (110) planes normal to the x- 
axis; the normal stress is shown for a crystal with (100) planes 
normal to the a;-axis as before (dashed line). 



creased monotonically with strain once plastic flow oc- 
curred, because of thermal softening. In both uniaxial 
and shear cases, there were regions around a few per- 
cent strain where the flow stress dipped because of plastic 
heating. A crystal would be susceptible to strain localiza- 
tion under these conditions. The shear stress increased 
with strain, but showed a clear dependence on strain rate 
and crystal orientation. The relative magnitude of flow 



FIG. 11: Effect of strain rate on shear stress in Ta crystals 
compressed uniaxially along [100] (solid) and [110] (dashed) 
directions. The dislocation density was held constant at 
W/fiia^ and there were no obstacles. 




Strain 



FIG. 12: Effect of strain rate on shear stress in Ta crystals 
sheared in [100] (solid) and [110] (dashed) planes. The dislo- 
cation density was held constant at 10//im^ and there were 
no obstacles. 



stress in [100] and [110] crystals changed with strain rate. 
(Figs [H and [H) 

It is interesting to compare these calculations with DD 
simulations and with experimental observations. One 
feature of the calculations is that initial flow stresses are 
in the range of around 0.7-3.5 GPa. The strain rate to 
achieve a given flow stress depends on the dislocation 
density. DD simulations show plastic yield at a similar 
stress, but with work-hardening at larger strains as dislo- 
cations tangle [s^j . Projectile impact experiments on cast 
and machined Ta samples millimeters thick have given 
low pressure flow stresses between 0.77 and 1.1 GPa fsof . 
Laser-induced shock experiments on rolled foils of Ta 25- 
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50 /im thick have given flow stresses of around 3 GPa [5l| . 
Thus the calculations presented here, though with a sim- 
plified model, are quite consistent with experiments on 
widely ranging time scales and microstructures. 



VI. DISCUSSION 

The theory and model described here include much of 
the physics behind plastic flow. The purpose of the model 
is to simulate the response of samples to dynamic loading, 
resolving the samples spatially so that each single crystal 
is represented typically by tens of numerical cells. This 
resolution treats the variation in loading experienced by 
different parts of each grain, and allows the variation in 
response to be calculated. However, in the continuum dy- 
namics implementation used so far, the evolution of the 
defect population is treated locally in each numerical cell. 
In reality, defects such as dislocations may move through- 
out the grain, and may also pass into adjacent grains. In 
grains with initially few pinning sites, dislocations may 
travel a long way. It would be an improvement to treat 
defect transport in the solution of the continuum dynam- 
ics equations. Similarly, grain boundaries may act as 
defect sources. Source terms are included in the present 
theory, and could be altered in cells in contact with grain 
boundaries, but transport would be needed to allow the 
defects to propagate further than the numerical cell size. 

The fraction of plastic work converted to heat, /p, also 
known as the Taylor- Quinney factor [52] , is an interesting 
and poorly-understood quantity. As discussed above, we 
have outlined a way to estimate 1 — /p from the potential 
energy of the dislocation population, but we have not yet 
investigated the properties of this approach in practice. 
The total dislocation core energy can be estimated from 
the dislocation density. It may be possible to quantify 
the energy associated with pins and other obstacles in a 
similar way. These contributions could then be included 
directly in the theory. The elastic strain energy associ- 
ated with the defects is likely to depend more sensitively 
on the configuration of the defects. It may be possible 
to estimate it from the dislocation density; the accuracy 
could be investigated by comparison with MD and DD 
simulations. 

The formulation chosen for defect-based flow provides 
an explicit, lattice-level mechanism for the Bauschinger 
effect through the dependence of the dislocation kinetics 
on the mass density, which differs between compression 
and release, though this is not the only contribution. De- 
fect evolution may be such that an initial plastic defor- 
mation reduces the flow stress on reverse loading, but this 
would depend on the relative behavior of the generation 
and annihilation terms and of the evolution of obstacle 
populations. If rotations are not excessive, it may well be 
that the initial deformation could increase the density of 
some dislocations without increasing the density of oth- 
ers which would act as obstacles, and then reverse load- 
ing would be accommodated in opposite-direction flow 



by the increased dislocation population. Contributions 
to the Bauschinger effect may also arise out of deforma- 
tion of a multi-grain microstructure. 

Although much of the plasticity theory presented here 
applies to twinning as well as slip, and although the 
disclinations separating twinned material can be viewed 
as a plane of dislocations, twinning is not fully repre- 
sented by this theory in its current form. At the grain 
level, twinning occurs by the conversion of material of 
one orientation to the twinned orientation. This re- 
orientation of the crystal axes means that the slip sys- 
tems change in real space, so any competing flow by slip 
should take this into account. If twinning is the domi- 
nant mode of deformation, and each grain is loaded such 
that only a single disclination mode is active in each, then 
re-orientation does not matter to the material response 
and the theory as presented should be able to treat the 
problem. 

We preferred to express the variation of elastic con- 
stants, the Peierls barrier, etc, in terms of the mass den- 
sity rather than the pressure. These quantities can be 
estimated from electronic structure calculations. Elec- 
tronic structure calculations using the local density ap- 
proximation for the exchange-correlation energy, which 
is almost universal, give results which typically have a 
systematic error o(l%) in mass density; this error should 
be corrected e.g. using the observed mass density at STP 
[3^ . A particular attraction in the use of electronic struc- 
ture calculations is that properties can be predicted for 
structural phases which are difficult to examine experi- 
mentally, such as those occurring at high pressures and 
temperatures induced in shock experiments. 

In most of the example simulations above, a strain was 
applied instantaneously and the relaxation of the stress 
state was predicted. In contrast, strength models are 
conventionally tested by applying a constant strain rate 
and predicting the flow stress. The plasticity theory pre- 
sented here is motivated by the need to analyze and in- 
terpret high pressure and strain rate experiments on sin- 
gle crystals. Typical experimental data are the velocity 
history at the surface of a planar crystal loaded with a 
square pressure history. For loading pressures slightly in 
excess of the Hugoniot elastic limit, the velocity history 
typically includes an elastic wave followed by a plastic 
wave. In single crystals, the rise time of the elastic wave 
is usually too rapid to measure [l^l, but the subsequent 
velocity history shows evidence of the time-dependence 
of plastic flow. Thus the strain rate during loading is very 
high and ill-determined; it is more appropriate in this in- 
stance to consider an instantaneous loading followed by 
relaxation at constant compression, as a preliminary to 
full continuum dynamics simulations in which the density 
may also change after the initial compression. 

At the high pressures and short time scales considered 
here, the bias effect of the high resolved stresses on the 
Peierls barrier can apparently allow shear stresses to re- 
lax much more rapidly. This effect is equivalent to a 
pressure-dependent reduction in the effective flow stress. 
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This is potentially important in the control of the seeding 

of implosion instabilities of microstructural origin in in- 
ertial confinement fusion: higher-pressure shocks may in- 
duce less spatial variation in material velocity caused by 
anisotropy in different grains. According to the simula- 
tions for Ta, it is also desirable to have a high dislocation 
density. 

The philosophy underlying physics-based models of 
material behavior is to have general models, parameters 
inherent to material of a given composition (e.g. Ta), 
together with characterization of the microstructure for 
each specific application. The characterization would in- 
clude the crystallographic texture, and defect and ob- 
stacle populations. The plasticity theory described here 
does not quite follow this ideal, as some microstructural 
information is at present implicit in the functional forms 
for the evolution of the defect populations. MD or DD 
studies are needed to elucidate the generation of defects 
and obstacles. 



VII. CONCLUSIONS 

We have developed a theory for elasticity and plas- 
tic flow intended to represent the response of crystals to 
deformation at high pressures and strain rates. The the- 
ory is needed because at pressures over 0(100 GPa) and 
on time scales shorter than 0(10 ns), elastic strains can 
exceed several percent, which changes the orientation of 
slip systems significantly with respect to the stress field. 
The motivation behind this work is to simulate the re- 
sponse of single crystal and polycrystallinc samples to dy- 
namic loading, by resolving the microstructure explicitly. 
The equations describing the response of the material are 
constructed to be convenient for integration forward in 
time. In this vein, we have preferred to consider the 
elastic strain rather than the stress as part of the local 
material state; this approach is more consistent with the 
conventional treatment of the scalar equation of state in 
continuum dynamics at high pressures. 

For a given total strain, the stress induced by the ma- 
terial was considered to originate from the elastic strain 
of the crystal lattice. Plastic flow was considered to op- 
erate at constant total strain, acting to relieve the elastic 
stress. With this approach, arbitrary sets of slip systems 
can be treated simultaneously, and in a spatially-resolved 
forward-time integration there is no need for explicit con- 
straints to ensure elastic compliance. If a crystal does not 
possess enough slip systems to relieve part of the stress, 
it remains unrelieved unless a grain elsewhere with a dif- 
ferent orientation can slip so as to relieve the total strain. 

We have proposed an underlying Arrhenius-based ki- 
netic term for the motion of dislocations and other de- 
fects which mediate plasticity. Defects do not 'know' 
which way to move: the applied stress was used as a bias 



on the Peierls barriers to motion, and the probability of 

hops in both directions (increasing as well as decreasing 
the elastic strain) was included. The barrier-lowering ef- 
fect of the stress was predicted to increase the plastic 
strain rate by hundreds of percent at low temperatures, 
and the reverse-hopping contribution was predicted to 
decrease the net strain rate by tens of percent at high 
temperatures. The barrier-lowering effect also allowed 
the ultimate theoretical strength of the material to be 
included in a self-consistent way; this was used in exam- 
ple simulations showing an abrupt increase in strain rate 
for suSiciently high applied stress. 

For the purposes of demonstration, we estimated the 
frequency factor in the Arrhenius rate for dislocation mo- 
tion from the Debye temperature, which is based on a 
simplified representation of the vibrational mode struc- 
ture in the crystal lattice, and depends on compression 
but not temperature. More rigorously, the frequency fac- 
tor could be calculated from the population of phonon 
modes which is able to promote a particular hop of a dis- 
location: a quantity which depends on temperature as 
well as compression. 

Although constructed with dislocations in mind, the 
treatment of defects is more general, and may be a suit- 
able framework for deformation by twinning. 

Ta was used as an example to illustrate the general 
behavior of the model, and to suggest ways to obtain pa- 
rameters for specific materials using more detailed theo- 
ries such as electronic structure, and experimental mea- 
surements of basic material properties. The calibration 
of parameters for any material is an involved process: de- 
tailed treatments of specific materials will be described 
separately. There is also much scope for future studies 
using molecular dynamics and dislocation dynamics to 
investigate the evolution of dislocations and pinning cen- 
ters. 
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